Minimum Energy Pathways via Quantum Monte Carlo 
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We perform quantum Monte Carlo (QMC) calculations to determine minimum energy pathways 
of simple chemical reactions, and compare the computed geometries and reaction barriers with those 
obtained with density functional theory (DFT) and quantum chemistry methods. We find that QMC 
performs in general significantly better than DFT, being also able to treat cases in which DFT is 
inaccurate or even unable to locate the transition state. Since the wave function form employed here 
is particularly simple and can be transferred to larger systems, we suggest that a QMC approach is 
both viable and useful for reactions difficult to address by DFT and system sizes too large for high 
level quantum chemistry methods. 
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Determining minimum energy pathways of reactions is 
of fundamental importance in scientific and technological 
applications. The knowledge of barrier heights (BH) is 
key to the prediction of catalytic properties of materials 
since it enables the use of transition state theory (TST) 
to determine reaction rates [IH3] ■ Locating efficiently the 
transition state on a potential energy surface (PES) is in 
fact a popular subject in computational physics and, to 
this aim, a variety of algorithms have been developed 
such as the shallowest ascent, synchronous transit, and 
nudged elastic band (NEB) approaches [4]. All these 
techniques ultimately rely on a method to determine the 
energy of a given atomic configuration and/or its deriva- 
tives. 

If we restrict ourself to quantum simulations, the most 
used approaches are density functional theory (DFT) 
or highly-correlated quantum chemical methods, that 
is, wave function post-Hartree-Fock techniques such as 
the coupled cluster single-double and perturbative triple 
approach [CCSD(T)], which is generally considered the 
"gold standard" in quantum chemistry. Many of these 
wave function methods are variational (though coupled 
cluster methods are not) and in principle offer a sys- 
tematic route to converge toward the exact energy, even 
though the increasing computational cost and the slow 
convergence severely limits this possibility. Their main 
drawback is that all these approaches implicitly or explic- 
itly rely on expanding the wave function in Slater deter- 
minants and, therefore, require large amount of computer 
memory and have a poor size scaling (iV 7 for CCSD(T), 
N being the number of electrons) , limiting their range of 
applicability to small systems. 

Consequently, for larger systems, DFT remains the 
method of choice due to its much more favorable com- 
putational cost (scaling from N 2 to iV 4 ). Even though 
continuous progress in the field has lead to the develop- 
ment of more precise and sophisticated DFT functionals, 
the situation is still far from satisfactory if one aims at 



high accuracy [SJ E] . For example, it is well known that 
popular functionals such as B3LYP [7HS] often lead to 
poor TS geometries and BHs [TOIEEI]. Since DFT meth- 
ods are not variational and do not offer a systematic way 
to improve their estimates, one has to resort to different 
approaches if better accuracy is needed. 

Alternatively, one can employ quantum Monte Carlo 
(QMC) methods, such as variational (VMC) and diffu- 
sion (DMC) Monte Carlo. These well-established ab- 
initio techniques take advantage of Monte Carlo integra- 
tion over the full Hilbert space. In particular, VMC is 
a stochastic way of calculating expectation values of a 
complex trial wave function, which can be variationally 
optimized. DMC provides instead, with a higher com- 
putational cost, a stochastic ground-state solution to the 
full Schrodinger equation, given a fixed nodal surface (us- 
ing the fixed-node approximation in order to avoid the 
notorious fermion sign problem). Because integrations 
are performed in the full Hilbert space, one can make 
use of non separable wave functions, with the explicit 
electron-electron correlation encoded in a so-called Jas- 
trow factor. This allows for noteworthy accuracy already 
using a simple and non memory-intensive single determi- 
nant Slater- Jastrow wave function. Although consider- 
ably more expensive than DFT methods (scaling as ./V 3 
with a much larger prefactor) , DMC generally offers bet- 
ter accuracy with respect to DFT. Furthermore, QMC 
methods possess a variational principle, which is a use- 
ful feature when one has to evaluate energy differences 
as in TST. From a computational point of view, QMC 
codes can be made to scale linearly with the number of 
cores and are not particularly memory demanding, mak- 
ing them suitable for today's massively parallel super- 
computers. Finally, QMC methods offer in principle the 
possibility to push the calculation up to a desired accu- 
racy by employing wave functions of increasing complex- 
ity (although, from a practical point of view, one is likely 
to adopt simple wave functions for intermediate-to-large 
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sized systems due to the increased computational cost of 
multi-determinant wave functions). 

In this article, we present nudged elastic band and 
climbing image calculations [HI [T3], where the geome- 
try optimization of all the NEB images is done fully at 
the QMC level. For some representative challenging re- 
actions from the NHTBH38/04 database QUE] and for 
a hydrogen transfer reaction [15] . we determine transi- 
tion state geometries and forward-reverse barrier heights 
within VMC and DMC, and compare our results against 
several current DFT functionals and other wave function 
methods. We demonstrate that VMC is able to locate re- 
action geometries with higher accuracy than DFT, while 
DMC outperforms DFT in evaluating BHs. Thus a sen- 
sible strategy, in terms of accuracy versus computational 
cost, is to calculate DMC BHs on VMC geometries. 

Methodology. The computation of the reactant 
and product geometries, the NEB calculations, and 
the saddle-point location through the climbing-image 
method are all optimization procedures over the total 
energy, although with different constraints. In particu- 
lar, we use here the Newton optimization method, based 
on the knowledge of the first and second derivatives of 
the energy. For the reactions reported here, these deriva- 
tives are from QMC calculations based on correlated sam- 
pling but, for larger systems, analytic calculation of QMC 
forces is possible and advisable in order to reduce compu- 
tational effort [TS]. Unless otherwise stated, we employ 
a single determinant Slater- Jastrow wave function: 

*(r,i?) = L» r (</> T ,r T , R)D l {(t) 1 ,r l , R)J{r, R), (1) 

where {r} and {R} denote the electronic and nuclear po- 
sitions, respectively. The Jastrow factor, J(r,R), explic- 
itly depends on the inter-particle coordinates, and in- 
cludes electron- electron, electron- nucleus, and electron- 
electron- nucleus correlation terms [17]. The Slater de- 
terminants, and D^, are constructed from the sets of 
molecular orbitals {</>^} and {4>^} for the up- and down- 
spin electrons, respectively. We employ scalar- relativistic 
energy-consistent Hartree-Fock pseudopotentials specifi- 
cally constructed for QMC calculations and expand the 
molecular orbitals on the corresponding cc-pVDZ basis 
set [THH5U]. These pseudopotentials have been exten- 
sively benchmarked, and their reliability has been re- 
cently supported by a DMC computation of atomiza- 
tion energies to near-chemical accuracy |21j . The pseu- 
dopotentials are treated beyond the locality approxima- 
tion [22]. 

Our choice of such minimal wavefunction and basis set 
is intentional since we want to maximize the scalability 
of our approach to systems larger than the ones consid- 
ered here. Our interest here is not to challenge quan- 
tum chemistry methods for small systems but, rather, 
to devise a strategy that has a more extended range of 
applicability while preserving a notable accuracy. While 
most DMC calculations found in literature use molecular 



orbitals computed with some other electronic structure 
method [23ti29], most often DFT, a key feature of our 
approach is that it is fully ab-initio since, at each iter- 
ation step in our geometric optimization, we perform a 
QMC optimization of all wave function parameters. This 
is done in order to guarantee consistency between the 
forces and the PES (see below) as well as to improve 
the results in terms of the absolute energy. The wave 
function optimization is performed at VMC level and de- 
tails of the procedure are described elsewhere [3D]. It 
is found that this optimization procedure only approxi- 
mately doubles the computer time needed to perform the 
calculations, while significantly lowering the expectation 
value of the energy. 

The QMC calculations are performed with a modified 
version of the CHAMP program [3T] . Since forces calcu- 
lated in QMC possess a statistical uncertainty, strictly 
speaking the optimization procedure via the Newton 
method never converges, so the final positions of the im- 
ages are obtained by averaging the geometries over sev- 
eral iterations after stationarity is reached. We use a time 
step of 0.01 a.u. in the DMC calculations. The DFT 
and Hartree-Fock (HF) all-electron calculations are per- 
formed with the GAMESS package [35], using Dunning- 
type Correlation Consistent triple-zeta basis sets, aug- 
mented with a set of diffuse function (aug-cc-pVTZ). 
No relativistic correction is included, differently from the 
QMC calculations which employ scalar-relativistic pseu- 
dopotentials. These corrections are however small for 
the light elements considered here and will not affect our 
results. 

Results. We select four challenging reactions from the 
NHTBH38/04 database [10] plus one hydrogen transfer 
reaction. As best estimates, we use the atomic geome- 
tries for the initial, final, and transition states reported 
in the database and computed through a quadratic con- 
figuration interaction with single and double excitations 
(QCISD) optimization. For these geometries, the BHs 
estimated with the Wl method (a complete basis set 
extrapolation over CCSD(T)) are also available. The 
reference data for H + OH — >• H2 + O are from ext- 
CAS+1+2+Q calculations [15] , 

We initially focus on the H + F 2 — >• HF + F reaction 
and collect the DFT and QMC data in Table [T] To mea- 
sure how much the geometries differ from the best es- 
timates, we calculate the RMS deviations of the inter- 
atomic distance among all atoms in the initial/final/TS 
configurations with respect to the corresponding best- 
estimate geometries. In the Table, a forward barrier (Vf) 
of zero means that the DFT functional finds no transi- 
tion state (i.e. the reactants are unstable) with the re- 
verse barrier being the reactant-product energy differ- 
ence. Many DFT functionals fail in finding any transi- 
tion state for this reaction, including the hybrid func- 
tional PBEO [33J, while B3LYP and M06 [31] retrieve 
a saddle point but with large deviations over the best 
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H + F 2 -> HF + F 


BE VMC VMC CAS DMC DMC CAS 


HF LSDA BLYP B3LYP PBE PBEO M06 


BHs (Kcal/mol) 


V) 2.27 6±1 Oil 2.2 ±0.5 Oil 
V r 106.18 126 ±1 113 ±1 114 ±1 108 ± 1 


0.0 0.0 1.4 0.7 0.0 0.0 2.7 
127.1 80.4 90.9 100.9 87.8 100.8 107.8 


RMS deviation (A) 


Reactants 0.008 0.007 
Products 0.002 0.008 
TS 0.028 0.013 


0.067 0.010 0.037 0.002 0.018 0.019 0.020 
0.016 0.018 0.020 0.009 0.017 0.004 0.001 
2.193 1.792 - - 0.216 



TABLE I. BHs and RMS of geometric deviations, H + F 2 — > HF + F reaction. The RMS is calculated over the 
deviation of the interatomic distances of all the atoms from the best estimate geometry. 



estimate TS geometry. In Table [TJ we also report the 
initial/final/TS geometries computed via VMC forces, 
where the uncertainty on the interatomic bonds due to 
the statistical noise on the forces is about 0.002 A. These 
geometries come from a fully VMC NEB and climbing 
image calculations. VMC is able to retrieve even at the 
single-determinant level the initial, the final and, espe- 
cially, the transition state geometry with much better 
accuracy than DFT. It is interesting to notice that, for 
geometrical data, VMC also performs better than the 
hybrid-meta M06 functional which is constructed to fit 
the BHs calculated on the best-estimate geometries from 
the database [34]. Clearly, this fitting procedure does 
not always guarantee that the actual TS retrieved by 
the functional is near the best-estimate one. Although 
the VMC geometries are rather accurate, the predicted 
reverse energy barrier (V r ) is markedly overestimated. 
Performing a DMC calculation on the VMC geometries 
retrieves better energy estimates, but the error on the re- 
verse barrier of about 8 Kcal/mol is still quite significant. 
We find that, to improve this energy barrier, it is not 
useful to reoptimize the geometries via DMC forces: The 
use of DMC forces does not alter significantly the geome- 
tries (within a statistical uncertainty on the TS geome- 
try of about 0.008 A) and, consequently, the estimated 
BHs either. In order to improve over these BH values, 
it is possible to take advantage of the variational princi- 
ple available in QMC and to resort to the use of multi- 
determinant al wave function. In this case, we employ a 
simple complete active space (CAS) wave function corre- 
lating three electrons in the three active orbitals relevant 
for the reaction and recompute the energy barriers over 
the VMC geometries obtained at the single-determinant 
level. The resulting VMC BHs are improved and the 
DMC values become very similar to the best estimates. 
Although the use of CAS wavefunctions is not readily 
applicable to larger systems, there exists scalable tech- 
niques to improve over single-determinant wave functions 
through the design of multi-determinantal size-extensive 
[35] or backflow wave functions [361 EZ] • 

We now return to the simple Slater- Jastrow wave func- 
tion, and consider the other reactions. In Figure [l] we 
compare the difference between the geometries obtained 
by VMC and various DFT functionals with respect to 
the reference ones obtained with QCISD, using the same 
measuring criterion as in Table [l| The functionals re- 
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FIG. 1. RMS of deviation of interatomic distances 
from the QCISD geometries (A). Distances are calcu- 
lated among all atoms involved in the reactions. 

ported in this figure are all hybrids or meta-hybrids and 
are the ones returning the smallest geometric/energy de- 
viation among the ones we tested on this set of reac- 
tions, which are the same ones listed in Table [T] The 
generalized-gradient-approximation functionals used in 
the overwhelming majority of DFT applications of CI- 
NEB, entail significantly larger deviations. For example, 
the deviation of PBE is two to four times larger than that 
of PBEO for TS geometries reported here, and the RMS 
BHs deviation of PBE on this set is more than twice 
than that of PBEO. For equilibrium geometries, VMC 
performs at the level of the hybrid functionals. For the 
TS, it typically returns more accurate geometries, often 
performing much better than DFTs. Notwithstanding 
that the M06 is actually fitted to reproduce BHs for the 
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FIG. 2. BHs deviation from best estimates calculated 
by QMC and DFT Methods (Kcal/mol). The dash- 
dotted line represents the DMC values obtained by employing 
the CAS wavefunction for the H + F 2 — s> HF + F instead of 
the single determinant one (full line), see text. 



NHTBH38/04 reactions, VMC still performs better than 
this functional in evaluating TS geometries. This may be 
again due to the M06 parametrization procedure, which 
does not guarantee the accuracy of the actual TS struc- 
ture calculated by the functional. We do not recalculate 
the geometries employing DMC forces because, from the 
test performed, DMC forces improves VMC TS geome- 
tries only slightly, as these are already notably accurate. 
Furthermore, DMC geometry corrections are barely re- 
flected in the calculation of the BHs, as these energies are 
second order in the deviation from the actual equilibrium 
points. For these reasons, we speculate that, for calculat- 
ing BHs, the calculations of geometries at VMC level and 
of BHs at DMC level is the most sensible choice regarding 
the trade-off between accuracy and computational cost. 

In Figure [2j we report the forward and reverse reac- 
tion barriers. While VMC (not shown) performs less 
accurately than hybrid functional, DMC calculated on 
the VMC geometries significantly improves the BHs of 
all reactions upon the VMC values, and performs at the 
level of the hybrid DFT approaches. In particular, our 
QMC procedure is more accurate than the hybrid func- 
tionals B3LYP and PBEO, while the only functional per- 
forming on average at the same level or slightly better 
is M06, which is actually fitted to reproduce exactly the 
NHTBH38/04 BHs. Therefore, we suggest that it is not 
guaranteed that, on a different set of reactions, the M06 



functional would still perform as well as QMC. Note that, 
if we exclude the problematic H + F2^ HF + F reaction 
(which can be fixed in DMC with the use of a simple CAS 
wave function as shown above), single-determinant DMC 
performs on average better than all DFT approaches, 
M06 included. The second most difficult reaction for 
DMC is H + N 2 -> OH + N 2 , which is also known to be 
strongly multi-determinantal in character, but where the 
DMC result is still better than DFT for both the forward 
and reverse barriers. 

Remarks on QMC Forces calculation. We compute 
here the VMC and DMC energy derivatives with cor- 
related sampling. While the procedure is relatively 
straightforward in VMC, the calculation of forces in 
DMC is more involved and approximations are generally 
used. For details about the DMC algorithm, we refer the 
reader to Ref. [35]. In both VMC and DMC, one should 
know the derivatives of the energy with respect to the 
wave function parameters and of the parameters with re- 
spect to the nuclear coordinates in order to calculate the 
forces correctly. In VMC, we avoid this problem by fully 
optimizing the wave function so that all the derivatives 
with respect to the parameters are zero. In DMC, the 
bias in the forces due to wave function optimization at 
VMC level has been found negligible in a few test cases. 

Conclusions. We investigated the possibility of per- 
forming full-QMC reconstruction of minimum energy 
pathways of chemical reactions. Full optimization of the 
wave function parameters is carried out during each iter- 
ation, so the employed technique is internally fully con- 
sistent. Geometric optimization of the minimum energy 
pathway and of the transition state is done at VMC level, 
with the obtained geometries being more accurate than 
DFT ones, especially for transition states. It also demon- 
strates the ability to correctly locate the TS in cases in 
which DFT fails in returning accurate geometries. At 
DMC level, the method displays very good performance 
in evaluating barrier reaction heights, comparing favor- 
ably even against hybrid functionals. Therefore, our ap- 
proach of calculating the geometries at the VMC level 
and the BHs at the DMC level is most effective as far as 
performance over computational cost is concerned: Cal- 
culating DMC geometries is very expensive and, in the 
tested cases, it does not improve significantly the esti- 
mates, while calculating DMC energies over VMC ge- 
ometries is much cheaper and still retrieves good results. 
Since the employed wave function is of the simple Slater- 
Jastrow type, this technique is scalable to larger sys- 
tems. Our results indicate that, for intermediate-sized 
system reactions where quantum chemistry methods are 
not computationally viable, the QMC approach may be 
the most accurate technique currently available. 

We acknowledge Dr. O. Valsson for assistance in the 
usage of CHAMP and GAMESS packages, and Prof. S. 
Baroni for helpful discussions. 
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